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ABSTRACT 



We analyze the impact of the choice rotation law on equilibrium sequences of relativistic differentially-rotating neutron stars in 
axisymmetry. The maximum allowed mass for each model is strongly affected by the distribution of angular velocity along the radial 
direction and by the consequent degree of differential rotation. In order to study the wide parameter space implied by the choice of 
rotation law, we introduce a functional form that generalizes the so called "j-const. law" adopted in all previous work. Using this new 
rotation law we reproduce the angular velocity profile of differentially-rotating remnants from the coalescence of binary neutron stars 
in various 3-dimensional dynamical simulations. We compute equilibrium sequences of differentially rotating stars with a polytropic 
equation of state starting from the spherically symmetric static case. By analyzing the sequences at constant ratio, r/|iy|, of rotational 
kinetic energy to gravitational binding energy, we find that the parameters that best describe the binary neutron star remnants cannot 
produce equilibrium configurations with values of T/\W\ that exceed 0.14, the criterion for the onset of the secular instability. 
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1. Introduction 

A neutron star (NS), during most of its life, is considered to 
be a stationary and rigidly rotating object, apart from a tiny 
lag between the rotation of the superfl uid component and that 
of the no rmal fluid and the crust (e.g. iBavm et al] (1 19691) and 
|Pines^Alpar(1985)). 

In fact a nascent neutron star which is born in a supernova 
event is likely to rotate differentially at first before its angular 
velocity distribution evolves toward a uniform rotation. There 
are several mechanisms that account for the redistribution of 
the angular momentum. O ne such mecha nism is the shear vis- 
cosity of neutron matter (ISawved Il989l) . An estimate of the 
timescale in which the viscosity damps a neutron star's internal 
shear motion has received a lot of attention in the literature (e.g. 
iCutler & LindbioinlJTgSTh '). For a typical neutron star we expect 
a timescale of 10 - 1 00 yrs. In the presence of magnetic fields, 
the magnetic br aking (|Spruitlll999 ) or th e magnetorotational in- 
stability (MRI) (lBalbus&Hawlevlll991h may drastically accel- 
erate the redistrib ution of angular momentum to the order of 1 s 
(iDuez et al ] l2006h . 

Dififerential rotation plays an important role in the beginning 
and at the end of the life of a NS. In ots early life, strong differen- 
tial rotation of a massive cor e in a supernov a may affect the col- 
lapse and bounce dynamics (iDimmelmeier et al. 2002: Ott et al. 
|2004|) . A newly-born neutron star with strong differential rota- 
tion may lose its stability in the course of the rotational (and 
thermal) evolution, and may collapse to a black hole or other 
type of compact star (quark star, hybrid star). This may give 
a unique neutrino signal and a strong gravitational wave emis- 
sion from the supernova explosion of massive stars (see, e.g.. 
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lOtt et al.l (l2007h and lTakiwak i & Kotake'feOlO*) for recent stud- 
ies on gravitational wave mechanism from stellar core collapse). 

An exciting way in which a neutron star's life can end 
is in the inspiral and consequent merger with a binary 
companion. Recent results of numerical 3-dime nsional sim- 
ulation s of neutron star b inary m ergers (e.g.. iRuffert et al.l 
fmd); 'Ruff^ert & JanO (l200lh : IShibata & Urvul (|2002l) ' 
Shibata & Taniguchj (120061): lA nderson et al.l OOOsI) : 
Baiotti et jlj (120081) : iGiacomazzo et al.i (l20ld) ) have shown 
that there are cases in which the remnant of the merger has 
a significantly high degree of differential rotation such that it 
can sustain a total mass co nsiderably la rger than that of a uni- 
formly rotating star (Baiott i~et al.ll2008l) . These hyper-massive 
neutron stars (HMNSs) remain stable over many dynamical 
timesc a les before collapsin g to a BH (e.g. Shibata & Urvfll 
( I2002h . iBaiotti et all (l2008h ). In order to study the effect of 
the centrifugal force in supporting the HMNS, most of these 
papers show the angular velocity profile of the merged object 
before the eventual collapse to a black hole. The typical profile 
extracted from simulations characteristically shows a plateau 
near the rotation axis and a certain distance from the axis shows 
a nearly power law behavior. The power index seems to differ 
from simulation to simulation, but generally does not agree 
with the so-called "j-const. law" dEriguchi & Mulleilil985h of 
rotation, which has been so far the only choice available to 
construct equilibrium models of differentially rotating stars. 

Models of rapidly rotating stars in general relativity have 
been studied since the 70s, when large numerical comput- 
ing facilities became available. The centrifugal deformation 
and general relativistic gravity make these investigations fully 
reliant on numerica l meth ods. Since the pioneering work of 
iButterworth & Ipse3 (Il976l) . these studies have included pro- 
gressively more sophisticated aspects, such as nuclear matter 
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EOS, degrees of differential rotation, magneti c fields and quite 
recently meridional flows (iBirkl etal.' '20ltf). We here name 
the fo llowing references of such studies: Butterworth & Ipse^ 
(119761) : Komatsu etal. (1989) ; Bonazzol a et all dl993p ; 
Stergioulas & Friedman, {1995) ; ,Baumgarte et al.l (|2000|) ; 



Ansorg et al. (2002). For more comprehensive collection of lit- 
eratures, see Stergioulas (2003"). 

However, it is surprising to note that all of these studies have 
used just a single kind of ansatz on the rotational angular veloc- 
ity profile (namely the "j-const. law", which include uniformly- 
rotating stars as a limiting case (Section l272l i). 

In this work we introduce a new parametrized functional 
form of the rotation profile that enables us to investigate a 
broader class of differentially rotating stars. In particular, our 
new rotation law reproduces different power laws in the outer en- 
velope of the angular velocity distribution (see Sectior lZ2l i with 
which it is possible to match the profile of the HMNS . To ana- 
lyze the impact of the law of rotation on the maximum allowed 
mass for the differentially-rotating neutron star we construct 
equilibrium sequences by imposing a fixed value of T/\W\ close 
to what is the classical limit for the onset of the secular insta- 
bility; rotating neutron stars are known to be destabilized by the 
effect of dissipative mechanisms like viscosity and gravitational 
radiation reaction (Chan drasekha r-Friedman-Schutz (CPS) i n- 
stabiHty; IChandrasekh^ (Il970t) : iFriedman & Schut j (Il978l) ). 
We show that the classical secular instability criterion is satis- 
fied only for a limited class of the rotation profiles. 

The paper is organized as follows. In Section |2] we briefly 
review the formalism of the equations and the numerical method 
used to solve them. Then we describe the new functional 
form of the rotation profile. In Section [3] we construct equilib- 
rium sequences of rotating stars using the new rotation profile. 
Discussions are given in the final section. 



transformed into convenient integral equations by using appro- 
priate Green's functions. 

The spatial linear velocity V of the flow with respect to an 
observer with zero angular momentum is given by 



y = /-Vsin0(n-w) 



(3) 



where Q - u'^ ju' is the angular frequency of the fluid measured 
in the asymptotic inertial frame. We next introduce the specific 
angular momentum^ j: 



J 



(4) 



Using these quantities, we can write the equilibrium equations 
for a stationary configuration as 

+ Bav - tV^aV + , ■' dA^ = (A = r, 6). (5) 



e + p 



1-^2 



1-jQ 



Assuming the fluid to be barotropic, the condition of integrabil- 
ity for Eq.(|5]l, can be written as 



l-jQ. 



(6) 



by introducing an arbitrary functional g, which is detailed in the 
next subsection. The first integral of motion for a stationary so- 
lution can be written as 



+ V + - ln(l - + I 8{n)dn = 0. 

P 2 



(7) 



The components of Einstein's equation and the first integral 
of hydrostatic e quatio n can be iteratively solved as described in 
lKomatsuetan(ll989h 



2. Formulation 

2.1. Equations for fluid and spacetime 

We construct configurations for an axisymmetric and stationary 
rotating perfect fluid in general relativity. The spacetime is as- 
sumed to be asymptotically flat and the flow is assumed to be 
circular (the velocity is only in t he azimuthal d irection). In this 
case we have two Killing vectors (Bardeen"l970') and can choose 
a coordinate system in such a way that the metric of the space- 
time is written as ('e.g.. lKomatsu et al.l (119891) X 

ds^ = -e^^dt^ + {dr^ + r^dO^) + e^^r^ sin^ e(dif - cjdtf, (1) 

where the spherical polar coordinates (r, 6, (p) are used. The met- 
ric potentials a,/3, v and oj are functions of r and only. Q The 
energy momentum tensor T"* of a perfect fluid is 

T"'' = (e + p)u"u'' + pg"'' (2) 

where e, p and u" are the total energy density, the pressure and 
the four velocity, respectively. 

The basic equations are: 1) Einstein's equation for the metric 
potentials Gat - ^^Tat, 2) rest mass conservation Vaipu") - 0, 
which is trivially satisfied under the present assumptions and 3) 
stress-energy conservation VhT"'' - Q. 

As described in Komat su et alJ (Il989l) . the components of 
Einstein's equation are cast into 4 equations for the potentials, 
three of which are elliptic partial differential equations. They are 

' We use geometrized units, c = G = 1. 



2.2. Rotation profile 

To exploit the first integral of the hydrostatic equation, 
we need to impose the integrability condition for Eq.(|6ll. 
Different choices of functional form g{Q.) would lead to various 
classes of differential-rotation profiles. All previous studies of 
differentially-rotating relativistic stellar models assume the sim- 
plest linear functional form, 

g(D.) = a2(Q, - Q), (8) 

where A is a constant with the dimensions of a length and is 
the angular velocity on the axis of rotation. 

This choice is termed as "j-const. law" since the Newtonian 
limit of the specific angular momentum / is that of "j -const. law" 
in Newtonian rotating stars (Eriguchi & Miiller 1985). 

We here introduce a more flexible form of the rotation profile 

as 

- Q") 

g{^) = , (9) 

1 - ^aHCi" - Q?) 

where a, Rq and Qc are constants. The coiTesponding specific 
angular momentum is 

= -^Q(Q'' - n^). (10) 



- This quantity is one definition of specifi c angular momentum i n 
general relativity. For different definitions, see lKozlowski et aTl il978h . 
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a 


-1 


-4/3 


-2 


-4 






^-3/2 






power law in envelope 


j-const. 


Keplerian 


v-const. 


HMNS 



Table 1. Choosing different a, it is possible to introduce differ- 
ent power law distributions of the angular velocity in the outer 
region of a star We here tabulate the Newtonian limits of par- 
ticular choices of a. "j-const." refers to constant specific angular 
momentum, "v-const." is for constant linear velocity. "HMNS" 
is the approximate power law of some of nu merical relativis- 
tic simulations with HMN S formation (e.g., IShibata & Urvul 
aiotti et al 



the sta r we adopt the "compactified" coordinate as in lCook et al.l 
(Il992h . to take into account the exact boundary condition at the 
radial infinity. 

For a given EOS, a functional of the rotation profile and max- 
imum density, the code requires information on "how rapidly" 
the star is rotating. To this end, it is more convenient to specify a 
parameter that measures the deformation induced by the rotation 
rather tha n the angular m omentum or the rotational frequency. 
We follow lKomatsuetal.1 ([1989) in fixing the ratio of the polar 
coordinate radius of the star to the equatorial radius. The physi- 
cal parameters of rotation as well as other physical quantities are 
then computed. 



The physical significance of the parameters a, Rq and D.^ in 
the profile is easily seen if we consider its Newtonian limit. Then 
the angular frequency is written as 



where R 
R » R(), 



2\s 



(11) 



rsinO. For R <s: Rq, we have Q ~ Qc, while for 



Uo 



(12) 



This means that the rotation profile has an inner plateau inside 
R ~ /?() and a power law envelope for R » Rq. The introduc- 
tion of the index a is an important feature of the new profile: as 
shown in Table 1, it is possible to reproduce different rotation 
laws by choosing different values for a. 

For some interesting cases (a e {-1, -2, -4)), it is possible 
to integrate analytically the expression for g{Q.) (see Appendix). 
It is important to note that for a = -1 the y-constant type of 
law is recovered - the two functional forms must agree in the 
limit of weak gravity. This particular case is marginally stable 
under Rayleigh's local stability criterion of axisymmetric insta- 
bility, which states that the specific angular momentum must not 
decrease outward in a stable star Values of a smaller than -1 
satisfy this condition (and hence are stable), a > also satis- 
fies this condition, but a star cannot be spun up rapidly enough 
to allow the appearance of the CFS instabiUty, because it easily 
reaches its mass-shedding limit. 

We now apply the new rotation law to construct equilibrium 
configurations of differentially rotating neutron stars and study 
their dependence on a. Each configuration is defined by five pa- 
rameters: the axis ratio rp/r^, the polytropic index A^, the maxi- 
mum density p„,ax and the two parameters of the rotation law, Rq 
and a. 



3. Results 

3. 1 . Numerical code 

We construct sequences of axisymmetric and differentially ro- 
tating objects using the code discussed in Komatsu et al. ( 1989). 
The code iteratively solves Einstein's equation and the first inte- 
gral of motion for a stationary fluid configuration. 

For all computation we used a grid with 600 points in the r- 
direction and 300 in the 6- direction. The radial grid is set up in 
such a way that the inner uniformly-distributed 300 grid points 
just cover the whole star (on the equator the 300th point from 
the origin always corresponds to the surface of the star). Outside 



3.2. Equation of State 

In the current study we use a polytropic equation of state, with 
the index = 1 to model neutron stars. 



(13) 



Typical neutron star masses and radi are recovered when we set 
K = 100 in geometrized units. We fix this parameter for all the 
sequences to have a crude model of the nuclear equation of state 
for neutron matter. In a subsequent paper we will investigate the 
dependence of the space of configurations on the EOS, using 
realistic cold and finite-temperature EOSs. 

3.3. Parameters 

To obtain an equilibrium model, we fix the parameters of the 
rotation profile, a and R{), then choose the maximum density 
Pmax and the coordinate axis ratio. In our code this set of param- 
eters uniquely determines an equilibrium configuration, whose 
physical characteristics such as mass, compactness and angu- 
lar momentum, angular frequency at the centre are computed 
once the solution is obtained. Rest mass and gravitatio nal mass 
are computed by using the standard formula (see e.g. iBardeenI 
(1970)). An important meas ure of rotation , the T/\W\ parameter, 
is introduced as in Komatsu et al ] (119891) . This is used both in 
Newtonian and general relativistic studies of rotating stars and 
defined as the ratio of rotational k i netic e nergy to gravitational 
binding energy (see lKomatsu et al](ll989h for details). It charac- 
terises the overall strength of rotation of the star Classical stud- 
ies suggest that when T/\W\ ~ 0.14, bar-shaped equilibria bifur- 
cate (into "Jacobi" or "Dedekind" type equilibrium sequences) 
and T/\W\ ~ 0.27 marks the onset of dynam ical instability of ax- 
isymmetric configurations (lTassoullll978l) . In general relativity 
numerical simulations point to a rather lower limit, T/\W\ ~ 0.25 
(.Manca et al.,.2007.) . 

3.4. Numerical equilibrium sequences 

For a given rotation profile and a maximum density, we start 
from a non-rotating model (the polar-to-equatorial axis ratio is 
equal to unity), and decrease the ratio to obtain stars with in- 
creasingly rapid rotation. We terminate the sequence if one of 
the following conditions is satisfied: 1) T/\W\ = 0.14. This is the 
classical criterion at which CFS instability or viscous instabil- 
ity sets in for a bar-shaped deformation (the actual critical point 
depends weakly o n both the equation of sta te and degree of dif- 
ferential rotation (Karino & Eriguchi 2002)). 2) Mass-shedding 
limit. The last condition is reached when the angular frequency 
of matter at the equatorial surface comes close to local Keplerian 
frequency and a further spin-up of a star leads to shedding of 



3 



Filippo Galeazzi, Shin'ichirou Yoshida, Yoshiharu Eriguchi: Differentially rotating neutron stars 



mass. Beyond this point it is not possible to construct equilib- 
rium configurations. 3) Topology change. When the degree of 
diff'erential rotation is large the maximum density may move 
away from the rotation axis. When the value of the axis ratio 
reaches zero, the star develops a hole on the rotation axis. The 
equilibrium sequence may be continued beyond this point, but 
we are not interested in this toroidal shape configuration as a 
model of a neutron star. 

Concerning the last point above, we should note that the 
structure of the parameter spac e is quite complex, as has been an- 
alyzed bv lAnsorg et alj (|2009|) . The toroidal and the spheroidal 
configurations form disjoint families of solutions separated by a 
mass shedding region. Increasing the degree of differential rota- 
tion by reducing the value of the /?() parameter, the two families 
join again in the parameter space. The area of the region separat- 
ing the two topologically different families depends on the p„,ax 
of the configuration, becoming smaller as the latter is increasing. 
Stated differently, starting from the Newtonian limit and going 
to the more general relativistic case, it becomes more and more 
complicated to find a value of Rq for which the solution can reach 
the toroidal family, in fact the separation increases between the 
two types of solutions. We have studied the parameter space of 
the solutions for the three different rotation laws analyzed here, 
and will collect the results along with the dependency hidden in 
the EOS in the followup paper. 

3.5. Stellar mass in parameter space 

To visualize our results, we plot gravitational mass of stellar 
models in the parameter space. Figure [1] shows the surfaces of 
equilibrium models embedded in the Rq - p„„ ^ - M/Msun space, 
where M is the gravitational mass in units of solar masses, 
for three different values of the parameter a: a - -1 (top 
panel), a = -2 (middle panel) and a = -4 (bottom panel). 
In order to have the same degree of differential rotation for all 
the models with a con stant R^), following the prescription of 
iBaumgarte et al.1 (l2000h . we define: 



^0 - Ro/Rc 



(14) 



where Rc is the circumferential radius. 

On the bases of these plots, we marked each parameter pair 
(Pmax, ^o) with different symbols to indicate how the equilibrium 
sequence ends. The triangles correspond to the parameters for 
which T/\W\ of the model reaches 0.14 (and the computations 
are stopped there). The dots correspond to the parameter pairs 
for which the sequence terminates at the mass-shedding limit. 
The solid lines are for the parameters for which the sequence of 
spheroids ends with a topological change. 

As mentioned in Section IT?! our sequences terminate at the 
three different criteria. Therefore one needs to be careful in read- 
ing these figures. It is important to stress that for the sequences 
terminating at T/\W\ = 0.14 and at the topological change, the 
value of mass plotted here is not the maximum mass of the se- 
quence. The solutions of the toroidal class do not possess a value 
for the maximum mass, instead this quantity can arbitrarily in- 
crease as the torus becomes thinner and thinner along the se- 
quence of equilibrium figures (see. Ansorg et al. (2003)). 

It is seen from Eq.(|9]i that the angular velocity profile be- 
comes uniform as — > oo. In this rigid rotation limit, the se- 
quence of = 1 polytropes is known to terminate with mass- 
shedding state before reaching T/\W\ - 0.14. The degree of 
differential rotation depends both on the power-law exponent a 
and on the radius Rq. We may regard it as weak when the pro- 
file is close to uniform rotation (i.e. /?() — » oo and/or a — » -oo). 



As in the case of uniform rotation, for stars with weak differen- 
tial rotation the centrifugal force at the equator on the surface 
reaches a value at which the star sheds mass before its T/\W\ 
value reaches 0.14. In this case all the sequences terminate at 
mass-shedding limit before reaching T/\W\ - 0.14 or the point 
of topology change. However a star with sufficiently strong dif- 
ferential rotation can store large rotational energy deep inside the 
star and allow the surface angular frequency to be smaller than 
that of the mass-shedding limit. Therefore a star with a = -1 
or -2 reaches the critical point T/\W\ = 0.14 before encoun- 
tering mass-shedding or topology change, provided that is 
small enough (that is a smaller core region of nearly uniform 
rotation). For these choices of a and sufficiently small ^o, we 
also see the appearance of topological change before the critical 
T l\W\ - 0.14 is reached (represented as the triangles in FiglT). 
Indeed, the sequences which we terminate at T l\W\ - 0.14 will 
eventually see the topology change if extended to faster rotation. 

An important observation here is that for a - -A none of the 
models reach T l\W\ - 0.14 before mass shedding occurs, and 
for a - -1, and -2 we need to set sufficiently small (i.e., 
enabling high degree of differential rotation) to have the critical 
T l\W\ before topology change or mass shedding occurs. 

4. Summary and Discussion 

We have introduced a new rotation profile to study equilibrium 
sequences of differentially rotating relativistic stars. Compared 
with the previous studies, which assume only one type of rota- 
tion profile, we are now able to investigate a broader class of 
rotating stars. 

As a first step towards systematic studies, we focus on the 
simplest neutron star model with a polytropic EOS with in- 
dex = 1. For the rotation profiles that allow analytic ex- 
pressions in the Eq. (|5]i, we computed sequences of spheroidal 
equilibria that start from non-rotating stars. Special attention 
is paid to the appearance of the critical point to the secular 
instability (CFS or viscous instability) to bar-shaped deforma- 
tion of the star. We see that the appearance of the critical point 
may occur for rather strong differential rotation in the case of 
a - -\ and -2 in Eq.®, which corresponds respectively to so- 
called "i;Const;[]_andJS;-coi^ rotation profiles in Newtonian 
stars (lEriguchi & Mulleii[T98l . Even in these cases, the critical 
point does not appear before the configuration changes from a 
spheroidal topology to a toroidal one, unless the parameter Rq in 
Eq.© is sufficiently small. 

When a - -A instead there seems no critical point of CFS 
instability before the equilibrium sequences reach their mass- 
shedding limit. This last case is relevant since it appears to 
mimic the approximate rotation profile of some of the quasi- 
stationary HMNSs seen in the numerical simulations of neutron 
star mergers. However these stars may be susceptible to the rel- 
atively new class of dynami cal instability (so-called "low 
instabihtY", e.g. IC entrella et al. (200lh: IShibata et al.' (2002^ 
Watts et al] (l200l : .Saiio & Yoshidal (l2006h : lOu & Tohhn3 



( 2006h : BaiottietaL ( 20081)), whose existence rehe s on th e 
strong shear flow ("Watts et al.1 (l2005h : |Cprvino et al.' (2010^). 
Indeed, numerical simulations of binary NSs do i ndica te that the 
HMNS deve l ops a bar deformation (e.g.'Shibat a & Urvul (l2002h : 
iBaiotti et al.l (l2008l) ) on top of the background star with approx- 
imate axisymmetry. It would be interesting to study the appear- 
ance of low T l\W\ instability by using our equilibrium stars with 
different rotation profiles since the survival time of the merger 
remnant could be an important diagnostic. The observation of 
the delay between detection of the gravitational wave signal 
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Fig. 1. Two-dimensional surfaces of equilibrium models of dif- 
ferentially rotating neutron stars embedded in the space - 
Pmax - M/Msun- The first panel is for a - -I, the second one 
is for ff = -2 and the last one is for a = -4. The lines at the 
bottom of each panel are drawn to clarify which of the three ter- 
minal conditions (see Section [34b of sequence applies. Triangles 
are for sequences reaching T/\W\ - 0.14, dots are for sequences 
terminating with mass shedding and thick solid lines are for the 
case when topology change occurs. 



form a binary neutron stars merger and the possible observation 
of the short gamma ray burst counterpart would then help to put 
constraints on the internal structure of a NS. 
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Appendix A: Integration of g{Q.) 

For some interesting cases, we have an analytic expression for 
the integral of g(Q) which appears in Eq.(|7]|. 

We define a parameter k := (RoQ.c)^ and normalize Q as 
^ := Q./Q.C- Then the integration of the functional giO.) is written 
as 



For a - -1, the right hand side reduces to 

1 



(A.l) 



-Arctan 



4-k 



4-k 

For a = -2, it reduces to 

\n(l+k(f - l))-2feln^ 
2{k-l) ■ 

Finally, for a = -4 we have 



(2^-1) 



■ln(l+^^(^-l)). (A.2) 



(A.3) 



-In^H- 



1 



In 



'-1 + Vl +4A:2-2fe^2^ 

^+1 + yf\TAi? + 2kf, 



(A.4) 
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